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ABSTRACT 

We study stability of gas accretion in Active Galactic Nuclei (AGN). Our grid based simulations 
cover a radial range from 0.1 to 200 pc, which may enable to link the galactic/cosmological simulations 
with small scale black hole accretion models within a few hundreds of Schwarschild radii. Here, 
as in previous studies by our group, we include gas radiative cooling as well as heating by a sub- 
Eddington X-ray source near the central supermassive black hole of 1O 8 M . Our theoretical estimates 
and simulations show that for the X-ray luminosity, Lx ~ 0.008 L B dd, the gas is thermally and 
convectivelly unstable within the computational domain. In the simulations, we observe that very 
tiny fluctuations in an initially smooth, spherically symmetric, accretion flow, grow first linearly and 
then non-linearly. Consequently, an initially one-phase flow relatively quickly transitions into a two- 
phase/cold-hot accretion flow. For Lx = 0.015 LiEdd or higher, the cold clouds continue to accrete 
but in some regions of the hot phase, the gas starts to move outward. For Lx < 0.015 Lsdd, the 
cold phase contribution to the total mass accretion rate only moderately dominates over the hot 
phase contribution. This result might have some consequences for cosmological simulations of the 
so-called AGN feedback problem. Our simulations confirm the previous results of Barai et al. (2012) 
who used smoothed particle hydrodynamic (SPH) simulations to tackle the same problem. However 
here, because we use a grid based code to solve equations in 1-D and 2-D, we are able to follow the 
gas dynamics at much higher spacial resolution and for longer time in comparison to the 3-D SPH 
simulations. One of new features revealed by our simulations is that the cold condensations in the 
accretion flow initially form long filaments, but at the later times, those filaments may break into 
smaller clouds advected outwards within the hot outflow. Therefore, these simulations may serve as 
an attractive model for the so-called Narrow Line Region in AGN. 

Subject headings: Galaxies: active, Galaxies: nuclei, Galaxies: Seyfert, Black hole physics, ISM: 
clouds 
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1. INTRODUCTION 

Physics within the central parsecs of a galaxy is domi- 
nated by the gravitational potential of a compact super- 
massive object. In_a classical theory of spherical accre- 
tion by iBondil (|1952D . Bondi radius Rb determines the 
zone of the gravitational influence of a central object and 
it is given by R B w 15O(Af BH /lO 8 M )(T oo /lO 5 K)- 1 pc, 
where Mbh is the central object mass, and is the 
temperature of the uniform surrounding medium. At 
radii smaller than the Bondi radius, Rb, the interstellar 
medium (ISM) or at least its part is expected to turn 
into an accretion flow. 

Physics of any part of a galaxy is complex. However 
near the Bondi radius, it is particularly so because there, 
several processes compete to dominate not only the dy- 
namical state of matter but also other states such as 
thermal and ionization. Therefore studies of the central 
parsec of a galaxy require incorporation of processes and 
their interactions that are typically considered separately 
in specialized areas of astrophysics, e.g., the black hole 
accretion, physics of ISM and of the galaxy formation and 
evolution. One of the main goals of studying the central 
region of a galaxy is to understand various possible con- 
nections between a supermassive black hole (SMBH) and 
its host galaxy. 

Electromagnetic radiation provides one of such con- 



nections. For example, the powerful radiation emitted 
by an AGN, as it propagates throughout the galaxy, 
can heat up and ionize the ISM. Subsequently, accre- 
tion could be slowed down, stopped or turned into an 
outflow if the ISM become unbound. Studies of heated 
accretion flows have a long history. Examples of early 
and key works includelOstriker et all (ll976luCowie et al.l 
(fl97l,lMathews fc Bregmad (119781) .[S tcllingwerf (1982). 
Bisnoyatvi-Kogan fe Blinnikovl (I1980D . iKrolik fe London! 
(1983?). and lBalbus fe Sokerl (|1989D . 

The accretion flows and their related outflows are 
very complex phenomena. It is likely that several 
processes are responsible for driving an outflow, i.e., 
not just the energy of the radiation, as mentioned 
above, but also for example, the momentum carried 
by the radiation. Therefore, our group explored com- 
bined effects of the radiation energy and momentum 
on the accretion flows and on producing outflows (e.g . 
Progall2007l; iProga et al.ll2008t Kurosawa fe Progall200S ; 
Kurosawa fc Progal l2009ri Kurosawa fc Progal 1200*9*3 : 
Kurosawa et al. These papers reported on results 



from simulations carried out using Eulerian finite differ- 
ence codes where effects of gas rotation and other com- 
plications such non-sp herical and non-azi muthal effects 
were included (see e.g. Uaniuk et al.ll2008l) . 

To identify the key processes determining the gas prop- 
erties (here, we are mainly concerned with thermal prop- 
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an accretion flow, in this paper we adopt a relatively sim- 
ple physical set up. Namely, the modeled system consists 
of a central SMBH of mass Mbh = 10 8 Af Q and a spher- 
ical shell of gas inflowing to the center. The simulations 
focus on regions between 0.1-200 pc from the central ob- 
ject, where the outer boundary is outside of Rb- The 
key difference compared to the Bondi problem is an as- 
sumption that the central accretion flow is a point-like 
X-ray source. The X-rays illuminate the accreting gas 
and the gas itself is allowed to cool radiativcly under op- 
tically thin conditions. To keep the problem as simple as 
possible, the radiation luminosity is kept fixed instead of 
being computed based on the actual accretion rate for an 
assum ed radiation efficiency (see also Kurosawa fc Progal 
I2009bl) . 

To model the presented problem, one needs to intro- 
duce extra terms into the energy equation to account for 
energy losses and gains. The physics of an optically thin 
gas that is radiatively heated and cooled, in particular, 
its thermal and dynam i cal sta bility has been analyzed in 
a great detail by iFieldl (|1965f ). Therefore, to study ther- 
mal properties of accretion flows or dynamical properties 
of the r mally unst able gas, it is worthwhile to combine 
iBondil (|1952t l and lFieldl (|1965t l theories. Notice, that our 
set-up is very similar that that used in the early works in 
the 70-ties and 80-ties that we mentioned above. Some 
kind of complexity and time variability in a heated ac- 
cretion flow is expected based on the 1-D results from 
the early work. 

A dynamical study of the introduced physical problem 
requires resolving many orders of magnitude of the ra- 
dial distance from the black hole. Our goal is not only to 
cover the largest radial span as possible but also to re- 
solve any small scale structure of the infalling gas. This is 
a challenging goal. To study the dynamics of gas in a rel- 
atively well controlled computer experiment , we use the 
Euler ian finite difference code ZEUS-MP ()Haves et al.l 
f2006h . 

We address systematically numerical requirements to 
adequately treat the problem of thermally unstable ac- 
cretion flows. We introduce an accurate heating-cooling 
scheme that incorporates all relevant physical processes 
of X-ray heating and radiative cooling. Low optical 
thickness is assumed, which decouples fluid and radia- 
tion evolution. We resolve three orders of magnitude in 
the radial range by using a logarithmic grid where the 
logarithm base is adjusted to the physical conditions. 
We solve hydrodynamical equations in 1- and 2- spatial 
dimensions (1-D and 2-D). We follow the flow dynamics 
for a long time scale in order to investigate the non-linear 
phase of gas evolution. Notice that most of the earlier 
work in the 70-ties and 80-ties, focused on linear analysis 
of stability, early stages of evolution of the solutions, and 
considered only 1-D cases. 

As useful our group's past studies are, we keep in 
mind that any result should be confirmed by using more 
than o ne techni que or approach. T herefore, iBarai et al.l 
(2012) (see also IBarai et al.l 120 1 If ) . began a parallel ef- 
fort to model accretion flows including the same physics 
but instead of performing simulations using a grid based 
code we used the s moothed partic le hydrodynamic (SPH) 
GADGET-3 code lSpringeJ ([20051). O verall, the 3-D SPH 
simulations presented by IBarai et al.l (|2012t ) showed that 
despite this very simple set up, accretion flows heated by 



even a relatively weak X-ray source (i.e., with the lumi- 
nosity around 1% of the Eddington luminosity) can un- 
dergo a complex time evolution and can have a very com- 
plex structure. However, the exact nature and robustness 
of these new 3-D results has not been fully established. 
IBarai et al.1 ()2012h mentioned some numerical issues, in 
particular, artificial viscosity and relatively poor spacial 
resolution in SPH, because of the usage of linear length 
scale (as opposed to logarithmic grid in ZEUS-MP), limit 
ability to perform a stability analysis where one wishes 
to introduce perturbations to an initially smooth, time 
independent solution with well controlled amplitude and 
spatial distribution (SPH simulations have intrinsic lim- 
itations in realizing a smooth flow). Therefore, the ro- 
bustness and stability of the solutions found in the SPH 
simulations are hard to access due to mixing of physical 
processes and numerical effects. Here, we aim at clari- 
fying the physics of these flows and measure the role of 
numerical effects in altering the effects of physical pro- 
cesses. 

Our ultimate goal is to provide insights that could help 
to interpret observations of AGN. We explore the condi- 
tions under which the two phase, hot and cold, medium 
near an AGN can form and exist. Such two phase accre- 
tion flows can be a hint to explain the modes of accretion 
observed in galactic nuclei, but also to explain the for- 
mation of broad and narrow lines which define AGN. We 
also measure the so-called covering and filling factors and 
other quantities in our simulations in order to relate the 
simulation to the origin of the broad and narrow line re- 
gions (BLRs and NLRs, respectively). The connection of 
this work to the galaxy evolution and cosmology is that 
we resolve lower spatial scales, and hence can probe what 
physical processes affect the accretion flow. In our mod- 
els, we can directly observe where the hot phase of accre- 
tion turns into a cold one or where an eventual outflow is 
laun ched. In most of the cu rrent simulations of galaxies 
(e.g. iDi Matteo et al.1120121 and references therein) these 
processes are assumed or modeled by simple, so called 
sub-resolution, approximations because, contrary to our 
simulations, the resolution is too low to capture the flow 
properties on adequately small scales. 

The article is organized as follows. In § [21 we present 
the basic equations describing the physical problem. In 
§[3j we show the details of the numerical set up. Results 
are in § [4] and in § O We summarize the results in § [6] 

2. BASIC EQUATIONS 
We solve equations of hydrodynamics: 

^ + PV.V = (1) 

Dv 

= -VP + pg (2) 

^( £ ) = -^V-v + p£ (3) 
Dt p 

where D/Dt is Lagrangian derivative and all other sym- 
bols have their usual meaning. To close the system of 
equations we adopt the P — (7 — l)e equation of state 
where 7 = 5/3. Here g is the gravitational acceleration 
near a point mass object in the center. The equation for 
the internal energy evolution has an additional term pC, 
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which accounts for gas heating and cooling by contin- 
uum X-ray radiation produced by an accretion flow near 
the central SMBH. The heating/cooling function con- 
tains four terms which are: (1) Compton heating/cooling 
(Gcompton), (2) heating and cooling due to photoioniza- 
tion and recombination (Gx), (3) free-free transitions 
cooling {Lb) and (4) cooling via line emission (L{) and it 
is given bv (IBlondinlli99l iProea et al.ll2000ft : 

pC = n 2 (G C om P ton + G X - L b - Li) [erg cm" 3 s^ 1 ] (4) 
where 
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where Tx is the radiative temperature of X-rays and T 
is the temperature of gas. We adopt a constant value 
T x = 1.16 x 10 8 K (E = 10 keV) at all times. The nu- 
merical constants in Equation |6] and [8] are taken from an 
analytical formula fit to the results from a ph otoioniza- 
tion code XSTAR (jKallman k Bautistal l2001). XSTAR 
calculates the ionization structure and cooling rates of 
a gas illuminated by X-ray radiation using atomic data. 
The photoionization parameter £ is defined as: 

fxLEdd fxLEddmpfi _, 

— = 5 — — ergs cm s J 

nr z pr z 
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where Fx is the radiation flux, n — p/(pm p ) is the num- 
ber density, and p, is a mean molecular weight. Given £ 
definition, notice that £ is a function of thermodynamic 
variables but also strongly depends on the distance from 
the SMBH. 

The luminosity of the central source Lx is expressed 
in units of the Eddington luminosities, fx = Lx/Lsdd- 
The reference Eddington luminosity for a supermassive 
black hole mass considered in this work is 
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[ergs s 
(10) 



3. METHOD AND INITIAL SETUP 



To solve Equation s. [T] [2 |3l_we use the numerical code 
ZEUS-MP (jHaves et alJl2006[ ). We modify the origi- 
nal version of the code in particular, we use a Newton- 
Raphson method to find roots of Equation [3] numerically 
at each time step. We have successfully tested the nu- 
merical method against an analytical model with heating 
and cooling. We describe the numerical code tests in the 
Appendix, showing the thermal instability (TI) develop- 
ment in the uniform medium. 

We solve equations in spherical-polar coordinates. Our 
computational domain extends in radius from 0.1 to 
200 pc. The useful reference unit is a radius of the 



innermost stable circular orbit of a central black hole: 
= &GMbh I c 2 ■ We assume the fiducial mass of the 
black hole M B h = 1O 8 M for which r* = 8.84 x 10 13 cm. 
The computational domain in these units ranges from 
n = 3484.27-, to r Q = 6.9683 x 10 6 r* (or r, = 6.610- 4 i? B 
and r a = 1.3Rb, where Rb — 152 pc). Since n is rela- 
tively large in comparison to the BH horizon we cannot 
model here the compact regions near the black hole where 
X-ray emission is produced. Instead we parameterize the 
X-ray luminosity using fx, so that Lx = fxLEdd- We 
solve equations for five values of /x=0.0005, 0.008, 0.01, 
0.015, 0.02 (these numbers correspond to models later 
labeled as A, B, C, D and E). 

As initial conditions for the A model (lowest luminos- 
ity), we use an adiabatic, semi-analytical solution from 
iBondil (|1952f h For higher luminosities the integration of 
equations starts from last data from a model with one 
level lower luminosity provided that the lower fx solu- 
tion is time-independent. The procedure is adopted in or- 
der to increase the luminosity in a gradual manner rather 
than sudden. Q Only for steady state solutions (with as- 
sumption that the mass accretion rate is constant from t\ 
to r.if ) the efficiency of conversion of gravitational energy 
into radiation r\ is related to fx as 



fx 
ra 



(11) 



where m is a mass accretion rate in Eddington units 

[Mem — LEdd/VrC 2 and rj r = 0.1 is a reference effi- 
ciency) and it is measured from the model data. In our 
steady state models, rhwl, therefore the energy conver- 
sion efficiency in these cases is approximately r/ — 0.1 fx- 
Our boundary conditions put constrains on a density 
at r a , it is set to be p Q — 10~ 23 gcm~ 3 . For other vari- 
ables we use an outflow type of boundary conditions 
at the inner and outer radial boundary. In 2-D mod- 
els our computational domain extends in £ (0,90°). 
At the symmetry axis and at the equator we use ap- 
propriate reflection boundary conditions. The numer- 
ical resolution used depends on the number of dimen- 
sions i.e.: in 1-D N r = 256,512, 1024,2048,4096; in 2-D 
(N r , N e ) = (256, 64), (512, 128), (1024, 256). The spacing 
of the radial grid is set as dri/dr i+i =1.023, 1.01, 1.0048, 
1.002, 1.0008 for A r =256, 512, 1024, 2048, and 4096, 
respectively. The number of grid points in the second 
dimension are chosen so that the linear size of the grid 
zone in all directions is similar (i.e., riA9j Ar,). 

4. RESULTS: 1-D MODELS 
4.1. 1-D Steady Solutions 

We begin with presenting the basic characteristics of 
1-D solutions. Table [T] shows a list of all our 1-D sim- 
ulations. Each simulation was performed until tf = 20 
Myr equivalent to 4.7 dynamical time scales at the outer 
boundary r Q = 200pc (t dyn = tff = v / rl/2GM B H = 

1 We also decouple Lx from M in order to avoid introducing 
additional parameters into the equations. While coupling these 
quantities not only a radiative efficiency of gravitational to radia- 
tive energy has to be assumed but one also needs to know how to 
calculate the mass accretion rate at the very compact region way 
below Ti = O.lpc. Another reason for decoupling Lx and M is that 
we are interested in caring out a stability analysis and perturb a 
steady state solutions with all model parameters fixed. 
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4.21 Myr). Only some of the numerical solutions settled 
down to a time- independent state at t/. We focus on ana- 
lyzing two representative solutions, that are steady-state 
at tf. 1D256C and 1D256D, with the X-ray luminosity of 
the former fx = 10~ 2 , and of the latter fx — 1.5 x 10~ 2 . 
Note that these solutions were obtained using the low- 
est resolution. We find these two solutions instructive in 
showing the thermal properties of the gas. 

Figure Q] presents the overall structure of model 
1D256 C and D (model C and D in the left and right 
column, respectively). Panels from top to bottom in Fig- 
ure [T] display: radial profiles of gas density, gas temper- 
ature overplotted with the Mach number (red line with 
the labels on the right hand side of the panels), the net 
heating/cooling rate plotted together with contribution 
from each physical process (see Equation!?]), the entropy 
S, and the bottom row shows gas temperature as a func- 
tion of £. In the bottom panels, the red line indicates 
the T-£ relation for radiative equilibrium (i.e. solving 
£(£, T) = for each T). The green line indicates a T — £ 
relation for a gas being adiabatically compressed due to 
the geometry of the spherical accretion (T oc £ 2 ), while 
the blue line for a constant pressure gas (T oc f 1 ). 

The 1D256C and D solutions differ mainly in the po- 
sition of the sonic point and in the fact that the model 
1D256D is strongly time dependent for a short period of 
time during initial evolution (see below). However, in 
most part the solution share several common properties. 
In particular, in both solutions, the gas is nearly in ra- 
diative equilibrium at large radii whereas, at small radii 
(below r sa 2 x 10 19 cm, where T > 2 x 10 6 K) they depart 
from the equilibrium quite significantly. In the inner and 
supersonic parts of the solutions, T scales with £ as if 
gas was under constant pressure. At large radii where 
the solutions are nearly in the radiative equilibrium, the 
net heating/cooling is not exactly zero. One can identify, 
four zones where either cooling or heating dominates. In 
the most inner regions where the gas is supersonic, adia- 
batic heating is very strong and the dominant radiative 
process is cooling by free-free emission. At the outer 
radii, the cooling in lines and heating by photoioniza- 
tion dominates. For models considered in this paper the 
Compton cooling is the least important. 

Inspecting the bottom panels in Figure [TJ one can sus- 
pect the gas is in the middle section of the computational 
domain to be thermally unstable because the slope of 
the T — £ relation (in the log-log scale) is larger than 1. 
Notice also that in both solutions the entropy is a non- 
monotonic function of radius. The regions where the en- 
tropy decreases with increasing radius correspond to the 
regions where there is net heating and the Schwarzschild 
criterion indicates convective instability at these radii. 
We therefore conclude that both solutions could be un- 
stable. We first check more formally the thermal stability 
of our solutions. 

4.2. Thermal Stability of Steady Accretion Flows 

The linear analysis of the growth of thermal modes un- 
der the radiative equilibrium con dition s (£(pg , Tp) = 0) 
has been examined in detail by iFieldl (|1965l) (see Ap- 
pendix for basic definitions). In Figure HJ in the top 
panels (left and right column correspond again to model 
1D256C and D), we show the radial profiles of various 




([ergs cm a -1 ] (l er g s cm s ~'] 

Fig. 1. — Structure of 1-D accretion flow in run 1D256C (left col- 
umn, fx = lxlO" 2 ) and 1D256D (right column, f x = 1.5xl0~ 2 ). 
Each panel is a snapshot taken at t=20 Myr. Panels from top 
to bottom show: density, temperature with Mach number (Mach 
number scale is on the right hand side), heating/cooling rates, and 
entropy S. The dashed vertical line in top panels marks the position 
of the sonic point. In panel with heating/cooling rates the black 
solid line is a net heating/cooling and color lines indicate partic- 
ular physical process included in the calculations: Compton heat- 
ing (red line), photoionization heating (green line), bremsstrahlung 
cooling (magenta-line), cooling through line emission (blue-line). 
The bottom panels display the gas temperature as a function of 
photoionization parameter; color lines indicate gas in radiative 
equilibrium (red), constant pressure conditions (blue) and free-fall 
compression (green). 

mode timescales. The timescalcs, t = 1/n, are calcu- 
lated using definitions I A2I and I A3! The growth timescale 
of short wavelength, isobaric condensations Tyj = —1/N P 
is positive (thermally unstable zone marked with the dot- 
ted line) in a limited radial range between about 10 and 
100 pc. The location of the thermally unstable zone de- 
pends on the central source luminosity, and it moves out- 
ward with increasing fx ■ The long wavelength, isochoric 
perturbations are damped, at all radii, on timescales of 
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Model ID 


fx 


N r 


*/ 
[Myr] 


(M)t 

[Mayr- 1 ] 


<X>r,t 


{fX,sc)t 


Max(r XiSC ) 


comment 


i D256 a 


5 x 10~ 4 


256 


20 


2.0 


3.9 


0.44 


0.49 




i dm 9 a 


5 x 10 — 4 


512 


20 


2 


6 1 


0.45 


0.51 




1 D1 094 A 


5 x 10 — 4 


1024 


20 


2 


6 8 


0.46 


0.52 




1 D904^ A 


5 x 10 — 4 


2048 


20 


2 


6 8 


0.46 


0.53 




1D4096A 


5 x 1CT 4 


4096 


20 


2.0 


6.8 


0.46 


0.53 


s 


1 F)9^fiR 


8 x 10 — 3 


256 


20 


1 8 


o 


1 


0.12 




inii OR 

_L 1 ' UlZi J_J 


8 x 10 — 3 


512 


20 


1.8 


o 


0.1 


0.13 




1 D1 024B 


8 x 10 — 3 


1024 


20 


1.8 


o 


0.1 


0.13 




1D2048B 


8 x 1CT 3 


2048 


20 


1.9 





0.1 


1.7 


s 


1D4096B 


8 x icr 3 


4096 


20 


1.95 


0.05 


0.13 


5.8 


11S 


1 D256C 


1 x 10~ 2 


256 


20 


1.7 


o 


0.09 


0.09 




1D512C 


1 x 10~ 2 


512 


20 


1.8 





0.09 


0.09 


s 


1D1024C 


i x icr 2 


1024 


20 


1.8 





0.11 


10.4 


s 




1 X 1U 


90A8 


90 
iU 




11 
U.Xi. 


n 19 
U.lo 


' -i 

y . o 


US 


1D256D 


1.5 x 10 -2 


256 


20 


1.5 


0.7 


0.2 


24 


s 


1D512D 


1.5 x icr 2 


512 


20 


1.5 


1.9 


0.23 


28 


IIS 


1D1024D 


1.5 x icr 2 


1024 


20 


2.1 


5.8 


0.55 


61 


IIS 


1D256E 


2 x 10~ 2 


256 


20 


1.23 


4.1 


0.2 


30 


IIS 



TABLE 1 

List of 1-D solutions. Columns from left to right are: model ID, fx - dimensionless luminosity, N r - numerical resolution 

IN RADIAL DIRECTION, tf - FINAL TIME OF THE SIMULATION, {M)t - TIME AVERAGED MASS ACCRETION RATE, \ - A RATIO OF COLD TO HOT 
MASS ACCRETION RATE AVERAGED OVER RADII LESS THAN 10 20 CM AND SIMULATION TIME, (COLD PHASE IS ANY GAS WITH TEMPERATURE 
T < 10 5 K), (T X> sc)e,t ~ ANGLE AND TIME AVERAGED OPTICAL DEPTH, MXx(t x ,sc) - MAXIMUM OPTICAL DEPTH RECORDED DURING THE 
SIMULATION, AND COMMENTS (S-REACHED A STEADY STATE SOLUTION AT tf , NS-NON STEADY STATE SOLUTION AT tf). 



t v — —1/N V (faster than TI development). The short 
wavelength nearly adiabatic, acoustic waves are damped 
as well, and r ac = —2/(N v — N p ). In FigureEJ the dashed 
line is the accretion timescale r acc = r/v. Within the 
thermally unstable zone, tti is short in comparison to 
T ace , in b oth models . 

IBalbusI (fl986lL IBalbus fc Sokerl (fl98l 
Math ews fc Bregmanl (|1978t ), ( and alsojKrolik fc London! 
1983) extended the analysis by [Field| (|1965t l to spherical 
systems with gravity, in more general case when initially 
the gas is not in the radiative equilibrium. Their ap- 
proximate solution gives the formula for linear evolution 
of the short wavelength, isobaric, radial perturbation 
as it moves with smooth back ground accretion flow 
(Equation 23 in IBalbusI 119861 or Equation 4.12 in 
Balb us fc Sokerj [T989). Since the two presented solutions 
are close to radiative equilibrium, the approximate 
for mula for the growth of a comoving perturbation given 
bv IBalbus fc Sokerl (|1989D reduces to 

where N p (r') is a locally computed growth rate of a short 
waveleng th, isobaric p erturbation as defined in the Ap- 
pendix or lFieldl ()1965l ). r' is radius where N p (r') < 0, and 
S s is an initial amplitude of a perturbation at some start- 
ing radius r — r s . Using Equationll2[ the isobaric pertur- 
bation amplification factors are 8/5 s ss 10 10 , 10 16 , 10 19 
and 10 33 , for models 1D256A, B, C, and D, respectively. 
Notice that these amplification factors are calculated 
for the asymptotic, maximum physically allowed growth 
rate, n = —N p , which might not be numerically resolved. 
To quantify the role of TI in our simulations we ought 



to address the following question. What is the mini- 
mum amplitude and wavelength of a perturbation in our 
computer models? The smallest amplitude variability is 
due to machine precision errors, ^machine ~ 10~ 15 (for 
a double float computations). The typical A of these 
numerical fluctuations are of the order of the numerical 
resolution, Arv The discretization of the computational 
domain affects the TI growth rates in our models in two 
ways: (1) the numerical grid refinement limits the size 
of the smallest fragmentation that can be captured; (2) 
the rate at which the condensation grows in the numer- 
ical simulations depends on number of points resolving 
a condensation. As shown in Appendix the perturba- 
tion of a given A has to be resolved by 20, or more, grid 
points. A wavelength Ao for which n = —0.9973 x N p is 
shown in Figure [3] together with Ar^ as a function of ra- 
dius for models with A^ r =256, 512, 1024, 2048, and, 4096 
grid points. In low resolution models we marginally re- 
solve Ao- We therefore expect the TI fragmentations to 
grow slower than theoretical estimates. Reduction of the 
growth rate due to these numerical effects even by a fac- 
tor of a few is enough to suppress variability because of 
the strong exponential dependence. 

Thermal mode evolution depends not only on the nu- 
merical effects but also other processes affecting the 
flow. Figure |4] shows the comparison of time scales of 
physical processes involved: the compression due to ge- 
ometry of the inflow and stretching due to accretion 
dynamics. We expect that any eventual condensation 
formed from the smooth background which leaves the 
thermally unstable zone, would accrete with supersonic 
background velocity. From the continuity equation, the 
co-moving density evolution is a balance of two terms 
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Fig. 2. — Left and right panels correspond, respectively to runs 
1D256 C and D. Top panels show the instability growth rates in 
comparison to the accretion time scale (r aC c = r/v, dashed line). 
The time scale for the short wavelength isobaric mode growth is 
displayed as the dotted line while the damping rate as solid line 
(tjv )• Other two lines show the long- wavelength isochoric mode 
damping rate tjv (heavy line) and the effective acoustic waves 
damping time scale Tjv„ (light line). Middle panel: The dashed 

line is r aC c and solid line is tgv = l/^BVi where u% v > is the 
Brunt — Vaisala oscillation frequency for a spherical system. Solid 
lines show the regions which are unstable convectivelly. The dotted 
line indicates region where 0Jg V < ant ^ oscillations are possible. 
Bottom panels: the derivative d In T/ci In £ as a function of radius 
is shown as a solid line. (dlnT/dln^) a[ j for an adiabatic inflow is 
marked as dotted line, and dashed line is the same derivative for 
radiative equilibrium conditions. Horizontal one indicates slope 
of 1. 




r [cm 



Fig. 3. — Grid spacing (red lines) in models with iV r =256, 512, 
1024, 2048, and 4096 points and Ao (black, dotted line) as a func- 
tion of radius in models 1D256 C (left panel) and D (right panel). 

(1/ p)(Dp/ Dt) = —2v/r — dv/dr, i.e., the compression 
and tidal stretching. The amplitude of condensation 



grows in regions where there is compression due to ge- 
ometry and decreases in regions where fluid undergoes 
acceleration - it stretches the perturbation. In the mod- 
els 1D256D and C interior of the TI zone, the evolution 
of the perturbation is dominated by compression because 
the compression time scale is the shortest. 
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Fig. 4. — Time scales in 1-D, stationary models 1D256 C (left 
panel) and 1D256 D (right panel): accretion time scale (race, 
dashed line), compression time scale (r c , solid line), tidal stretch- 
ing time scale (t s , dotted-dashed line), and condensation growth 
time scale (tjv , dotted line). 



4.3. Convective Stability of Steady Accretion Flows 

In this subsection, we examine in more detail convec- 
tive stability of our solutions. In Figure [U (middle pan- 
els), we compare the accretion time scale r acc and the 
Brunt — Vaisala time scale tbv = associated with 

LOBV 



the development of convection. The frequency ujbv is 
defined as u| v = The convectivelly un- 

stable regions are marked as solid lines (u| v > 0). The 
convectivelly unstable zones overlap with the thermally 
unstable zones. Since r acc <C tbv convective motions 
might not develop, at least at the linear stage of the de- 
velopment of TI. 

In the bottom panels of Figure [2] we show the log- 
arithmic derivatives of dlnT/dln^ that could be used 
to graphically assess the stability of the flow. This can 
be done by comparing the derivatives (the slopes of the 
InT— In £ relation) for three cases: model data (solid line 
where T and £ are taken directly from the simulations), 
purely adiabatic inflow (dotted line, assuming that the 
velocity profile is same as in the numerical solution), and 
radiative equilibrium conditions (dashed line). In partic- 
ular, the regions where the solid line is above the red 
line correspond to the potentially TI zones. The regions 
where the dotted line is below the solid line correspond 
with the zone where the flow is potentially convectivelly 
unstable. The conclusion regarding the flow stability is 
consistent with the conclusion reached by analyzing the 
time scales shown in the top and middle panels of Fig. 
2. 

4.4. Other Physical Consequences of Radiative Heating 
and Cooling - obscuration effects 

The growth of the thermal instability leads to the de- 
velopment of a dense cold clouds (shells in 1-D models; 
e.g. variable phase in model 1D256D). The enhanced ab- 
sorption in the dense condensations may make them op- 
tically thick. Here we check if the time-dependent mod- 
els are self-consistent with our optically thin assumption. 
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Fig. 5. — Fraction of central illuminating source radiative energy 
intrinsically absorbed (upper panels) and emitted (bottom panels) 
by gas per second as a function of time. We show the steady state 
solutions 1D256C and 1D256D in left and right panels, respec- 
tively. Solid lines show the net absorption/emission and dashed 
lines indicate the intrinsic absorption and emission. 

Figure [5] shows the amount of energy absorbed and emit- 
ted by the gas (heating and cooling rates integrated over 
a volume at each time moment) in comparison to the 
luminosity of the central source in models 1D256C and 
D. Solid lines show the net rate of the energy exchange 
between radiation and matter (cooling function C inte- 
grated over the simulation volume) while dashed lines in- 
dicate the intrinsic absorption and emission (heating and 
cooling terms used in C are integrated independently). 
The net heating-cooling rate is mostly much lower than 
unity reflecting the fact that in the steady state the gas is 
nearly in a radiative equilibrium. During a variable phase 
(part of model 1-D256D) the energy absorbed by the ac- 
cretion flow (black solid line) becomes comparable to the 
X-ray luminosity of the central black hole. During this 
variable phase the optical thickness of accreting shells 
can increase up to tx,sc ~ 20 where the majority contri- 
bution to opacity is due to photoionization absorption. 
The average optical thickness increases in models with 
higher resolution indicating that the flow is more vari- 
able and condensations are denser. This increase in opti- 
cal depth is related to shells condensating much faster in 
runs with higher resolutions. The dense condensations 
falling towards the center could reduce the radiation flux 
in the accretion flow at larger distances. It is beyond the 
scope of the present paper to investigate the dependence 
of the flow dynamics on the optical thickness effects and 
we leave it to the future study. 

Significant X-ray absorption is related also to transfer 
of momentum from radiation field to the gas. To estimate 
the importance of the momentum exchange between ra- 
diation and matter, one can compute a relative radiation 
force: 

f force = — fx (13) 

(JTH 

where ax is the energy averaged X-ray cross-section. The 



momentum transfer is significant when f force > 1- Us- 
ing our expression for the heating function due to X- 
ray photoionization cjx/vth = Hx/n/Fx — 2.85 x 

lfJ 4£-3/4 T -l/2 ( gee § Eyen for a dense cold shell 

f force is at most 0.1 (in case when T max ss 60). Therefore 
radiation force is not likely to directly launch an outflow. 
However, the situation may change when optical effects 
are taken into account. 

4.5. M Evolution 

We end our presentation of 1-D results with a few com- 
ments on the time evolution of the mass accretion rate, 
M. Figure [6] displays M vs time measured for all of our 
1-D models. One can divide the solutions into two sub- 
categories: steady and unsteady where M varies from 
small fluctuations to large changes. For a given fx, the 
time behavior of the solution depends on the resolution, 
due to effects described above. In the variable models 
a fraction of the accretion proceeds in a form of a cold 
phase defined as all gas with T < 10 5 K. Column 6 in Ta- 
ble Q] shows the ratio of cold to hot mass accretion rates, 
X, computed by averaging M's over the radius and sim- 
ulation time. We average M's over r < 100 pc because 
the cadence of our data dumps is comparable to the dy- 
namical time scale at 100 pc. The larger the luminosity 
the more matter is accreted via the cold phase. However, 
the maximum value of \ is of the order of a few, so the 
dominance of the cold gas is not too strong. 
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Fig. 6. — Mass accretion rate in 1-D solutions for fx = 
0.0005,0.08,0.01, and 0.015. Different colors shows M for various 
number of grid points: 7V r =256 (blue), 512 (red), 1024 (green), 
2048 (magenta), 4096 (black). 



5. RESULTS: 2-D MODELS 
5.1. Seeding the TI 

To investigate the growth of instabilities in 2-D, we 
solve eqs. Q] HI El for the same parameters Mbh and fx 
as in § 21 but on a 2-D, axisymmetric grid with 9 angle 
changing from to 90 deg. We use three sets of numerical 
resolutions described in § [3] To set initial conditions in 
axisymmetric models, we copy the solutions found in 1-D 
models onto the 2-D grid. In case of time-independent 1- 
D models, our starting point, is the data from t = tf. We 
checked, for example, the runs 2D256x64 A, B, C, and D 
(which are 2-D version of models 1D256A, B, C, and, D) 
are time-independent at all times as expected. In case 
of higher resolution models, for which the 1-D, steady 
state models do not exist we adopt a quasi-stationary 
data from the 1-D run early evolution (at t of a frac- 
tion of a Myr), during which the flow is already relaxed 
from its initial conditions but the TI fluctuations are not 
yet developed. Models which are time varying in 1-D, 
in 2-D develop dynamically evolving spherical shells, as 
expected, also indicating that our numerical code keeps 
a symmetry in higher dimensions. 

To break the symmetry in 2-D models, we perturb the 
smooth solutions adopted as initial conditions. The per- 
turbation of a smooth flow is seeded everywhere and 
has a small amplitude randomly chosen from a uni- 
form distribution. The new density at each point is 
p = po(l + Amp* rand) , where rand is a random number 
rand £ (—1,1) and maximum amplitude Amp = 10~ 3 . 
To seed the isobaric, divergence free fluctuations other 
(than p) hydrodynamical variables are left unchanged. 
The amplitude magnitude Amp is chosen to be much 
higher in comparison to e mac hine, in order to investigate 
the development and evolution of strongly non-linear TI 
on relatively short time scales, starting directly from a 
linear regime. The list of all perturbed, 2-D models is 
given in Table [2j 

5.2. Formation of Clouds, Filaments, & Rising Bubbles 

For luminosities Lx < 0.015 LEdd, the 2-D models 
show similar properties to the 1-D models. The gas is 
thermally and convectivelly unstable within the compu- 
tational domain, and we observe that very tiny fluctua- 
tions in an initially smooth, spherically symmetric, accre- 
tion flow, grow first linearly and then non-linearly. Since 
the symmetry is broken the cold phase of accretion forms 
many small clouds. For Lx = 0.015 LEdd or higher, the 
cold clouds continue to accrete but in some regions a hot 
phase of the gas starts to move outward. 

In Figure [71 we show three snapshots of representa- 
tive 2-D model 2D512xl28D at various times (t = 3, 6 
and 11.8 Myrs). This model has the best resolution and 
the highest luminosity for which we are able to start the 
evolution from nearly steady state conditions. Columns 
from left to right show density, temperature, and total 
gas velocity overplotted with the arrows indicating the 
direction of flow. Initially (at t=3 Myr) the smooth 
accretion flow fragments into many clouds, which are 
randomly distributed in space. The cooler, denser re- 
gions are embedded in a warm background medium. The 
colder clouds are stretched in the radial direction and 
they have varying sizes. This initial phase of the evolu- 
tion is common for all models in Table [2j 



The phase where many cold clouds accrete along with 
the warm background inflow is transient. At a later stage 
(t=6 Myr, middle panels), model 2D512xl28D shows a 
systematic outflow in form of rising, hot bubbles. The 
outflow is caused by the pressure imbalance between the 
cold and hot matter and buoyancy forces. The hot bub- 
bles expand at speeds of a few hundreds km/s. Despite 
of the outflow, the accretion is still possible. During 
the rising bubble phase, the smaller clouds merge and 
sink towards the inner boundary as streams/filaments. 
However, even this phase is relatively short-lived. Bot- 
tom panels in Figure [7j show the later phase of evolution 
when some of the filaments occasionally break into many 
clouds (this process takes place between 10 and 50 pc). 
These 'second generation' clouds occasionally flow out 
together with a hot bubble. Along the X-axis, we see an 
inflow of a dense filament. 

To quantify the properties of clumpy accretion flow, 
we measure the volume filling factor of a cold gas f vo i, 
defined as: 

fvol = V c l ou d/Vt t (14) 

where V c ioud is the volume occupied by gas of T < 10 5 
K, and Vtot is the total volume of the computational do- 
main. In model 2D512xl28D, the time- averaged fvoi is 
(fvol) = 3 x 10 -3 . The time evolution of f vo i within 
60 pc is shown in Figure |8] (black, solid line). The fvoi 
is variable and at the moment of the outflow formation, 
fvoi suddenly decreases by a factor of about 4. For com- 
parison fvoi calculated during run 2D256x64D is also 
shown (blue, dashed line). Run 2D256x64D has the same 
physical parameters as 2D512xl28D, however, no outflow 
forms. In the latter case, f vo i is less variable and larger. 
In Tabled we gather the time averaged (fvoi) for all 2-D 
solutions. Measuring fvoi allows to quantify whether the 
perturbed accretion flow returns to its original, smooth 
state. We find that this happens when the (fvol) ~ 10~ 5 
or smaller (models 2D256x64A, B, and C). 

5.3. M Evolution 

Figured] presents M through the inner boundary mea- 
sured as a function of time. In most cases (except 
model 2D256x64A), M becomes stochastic instantly with 
spikes corresponding to the accretion of colder but denser 
clouds similar to those in Figure [Jj (upper panels) . Sim- 
ilar to 1-D models, one can divide the solution into two 
types: steady and unsteady state. In the latter, M fluc- 
tuates on various levels depending on fx- 

Table [2] lists several characteristics of our 2-D simu- 
lations, for example, the ratio (x)t averaged over time. 
In 2-D models this variable is smaller in comparison to 
1-D due to geometry of the clouds. The maximum value 
of (x)t is less than unity. This indicates that multi- 
dimensional effects (specifically development of convec- 
tion) promote hot phases accretions. We plan to investi- 
gate this issue in future by earring out 3-D simulations. 

We find that a large scale outflow forms only in run 
2D512xl28D. But even in this case, the M is not signifi- 
cantly affected by the outflow. Figure fTOl shows the mass 
outflow rate (dashed line), inflow rate (dotted lines) and 
total mass flow rate (solid line) as a function of radius. 
The same types of lines show M for various times of the 
simulation (t=3, 6 and 11.8 Myr, green, blue and black 
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Model ID 


fx 


N r 


Ng 


tc 
[Myrl 

L J J 


(M) 
[Moyr- 1 ] 


\J V ol I t 


(y)r i 


(ty o.- ) a + 
\ -A. ,sc / tf , t 


Maxfrir B r-) 


final state 


2D256x64A 


5 x 10~ 4 


256 


64 


20 


2.0 


1 x 10~ 6 





0.45 


0.46 


smooth 


2D256x64B 


8 x 10~ 3 


256 


64 


15.4 


2.04 


1 x 10~ 6 


10~ e 


0.1 


0.99 


smooth 


2D512xl28B 


8 x 1CT 3 


512 


128 


20 


1.95 


1 x 10" 4 


0.02 


0.11 


4.7 


clouds 


2D1024x256B 


8 x 10~ 3 


1024 


256 


1.83 


1.84 


1.5 x 10~ 3 


0.39 


0.19 


24. 


clouds 


2D256x64C 


1 x 10" 2 


256 


64 


11.8 


1.94 


7 x 10-5 


0.15 


0.11 


2.5 


smooth 


2D512xl28C 


1 x 10~ 2 


512 


128 


20 


1.88 


5 x 10~ 4 


0.09 


0.13 


17.3 


clouds 


2D1024x256C 


1 x 10~ 2 


1024 


256 


1.12 


1.95 


4 x 10" 3 


0.5 


0.24 


70 


clouds 


2D256x64D 


1.5 x icr 2 


256 


64 


12 


1.57 


5 x 10" 3 


0.3 


0.14 


37.8 


clouds 


2D512xl28D 


1.5 x 10" 2 


512 


128 


11 


1.6 


3 x 10" 3 


0.43 


0.12 


13.2 


outflow, filaments 



TABLE 2 

List of 2-D solutions. Columns from left to right are: model ID, fx - dimensionless luminosity, N r ,Ng - numerical 

RESOLUTION IN RADIAL AND TANGENTIAL DIRECTION, tf - FINAL TIME OF THE SIMULATION, (M)t - TIME AVERAGED MASS ACCRETION RATE, 
(fvol)t ~ TIME AVERAGED VOLUME FILLING FACTOR, \ - A RATIO OF COLD TO HOT MASS ACCRETION RATE AVERAGED OVER RADII LESS THAN 

10 20 CM AND SIMULATION TIME, (COLD PHASE IS ANY GAS WITH TEMPERATURE T < 10 5 K), (rx,sc)e,t " ANGLE AND TIME AVERAGED 
OPTICAL DEPTH, MAX(t x , B c) " MAXIMUM OPTICAL DEPTH RECORDED DURING THE SIMULATION, AND COMMENTS. In THE BEGINNING OF ALL 
RUNS THE ACCRETION FLOW BECOMES UNSTEADY DUE TO RANDOMLY SEEDED PERTURBATIONS, AT tf THE FLOW EITHER RETURNS TO THE 
ORIGINAL, UNPERTURBED, SMOOTH STATE (GLOBALLY STABLE SOLUTIONS), OR REMAINS UNSTEADY WITH COEXISTING, TWO PHASE MEDIUM 
(COLD CLOUDS EMBEDDED IN A WARM GAS). ONLY IN ONE MODEL A LARGE SCALE OUTFLOW FORMS (MODEL 2D512xl28D). 



lines, respectively) and they are averaged over 9 angle. 
The rising bubble originates at about 10 pc in this case. 
We anticipate that large scale outflow are common and 
significant for high luminosity cases (i.e., for fx > 0.02). 

5.4. Obscuration Effects 

Here we again check whether the cloud opacity might 
affect our results. The averaged optical thickness of the 
filaments and clouds is similar (see Table [2] columns 8 
and 9). In Figure [TT] we show how much energy is ab- 
sorbed and emitted in run 2D512xl28D during the evo- 
lution. The figure shows the intrinsic absorption and 
emission integrated over the entire computational do- 
main. We next calculate {tx,cs) (optical thickness due 
to absorption, averaged over angles and times) and max- 
imum value of tx,cs that occurred during the evolution. 
In case of the largest optical depth of r « 70 (in run 
2D1024x256C) the radiation force coefficient from Equa- 
tion [13J f force ~ 0.2 which, as in 1-D models, is small 
but might not be negligible. Therefore, we are planning 
to explore the effects of optical depth in a follow-up pa- 
per. 

6. SUMMARY AND DISCUSSIONS 

In this work we show the evolution of thermal instabil- 
ities in gas accreting onto a supermassive black hole in 
an AGN. A simplified assumptions made in this work, 
in particular constant X-ray luminosity emitted near 
the central SMBH regardless of the M, allows to fol- 
low the development of TI from the linear to strongly 
non-linear and dynamical stage up to luminosities of 
L a 1.5 x 10- 2 L Edd . 

In our 1-D models the TI is seeded by numerical er- 
rors which might be non-isobaric and are initially under- 
resolved. In the initial phase, the TI growth rate is 
smaller than predicted by theory. The rate is affected 
by grid resolution which leads to the formation of cold 
clouds of various sizes and density contrasts. This is 
reflected in the mass accretion rate fluctuating at dif- 
ferent amplitude and rate for the same physical condi- 
tions but different resolutions. One cannot avoid dealing 
with these numerical difficulties in the numerical mod- 
els. Nevertheless, we find the under-resolved, 1-D models 



very useful in quick checking where the thermally unsta- 
ble zone exists and what type of fluctuation could cause 
the smooth to turn into a two-phase medium. For given 
physical conditions, Figure [3] shows the wavelength Ao 
and Equation [12] gives the amplitude of an isobaric per- 
turbation required to break the smooth flow into a two- 
phase, time-dependent model. 

In 2-D models, although the models depend on the res- 
olution effects same way as in 1-D setup, we can observe 
an outflow formation. The convectivelly unstable gas 
buoyantly rises and, as found in this work, controls the 
later evolution of the two-phase medium and mass accre- 
tion rate. Given a simple set up with minimum number 
of processes included, our models display the three major 
features needed to explain some of the AGN observations: 
cold inflow, hot outflow and cold, dense clouds which oc- 
casionally escape, advected with the hot wind. We show 
that an accretion flow at late, non-linear stages, thus 
most relevant to observations, are dominated by buoy- 
ancy instability not TI. This suggests that the numerical 
resolution might not have to be as high as that needed 
to capture the small scale TI modes and it is sufficient to 
capture significantly larger and slower buoyancy modes. 
We plan to check the consistency of the models with the 
observations by calculating the synthetic spectra, includ- 
ing emission and absorption lines base d on our simulatio n 
following an approach like the one in ISim et al.1 (12011) . 
Here, we only briefly comment on the main outflow prop- 
erties and compare them some observations of outflows 
in Seyfert galaxies. 

Space Telescope Imaging Spectrograph (STIS) on 
board the Hubble Space Telescope allows us to map the 
kinematics of the Narrow Line Regions in some nearby 



Seyfert Ga laxies (e.g . for NGC 4151 iDas et al.l 
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NGC 106 8: IDas et alBiioi M rk 3. ICrenshaw etal 
Mrk 573, iFischer et al.ll2010t and Mrk 78. iFischer et al.1 
2011). Position-dependent spectra in [O III] A 5007 and 
H a A 6563, and the measurements of the outflow velocity 
profiles show the following general trend: the outflow has 
a conical geometry and the [O III] emitting gas acceler- 
ates linearly up to some radius and then decelerates. The 
velocities typically reach up to about 1000 km/s and a 
turnover radius is on one hundred to a few hundred par- 
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Fig. 7. — Density, temperature, velocity magnitude with velocity direction (left, middle, and right column, respectively) in model 
2D512xl28D at t=3, 6 and 12 Myr (upper, middle, and bottom rows, respectively). Each panel show inner parts of the flow at r < 100 pc. 



sec scales. 

To compare our results with the observations, Fig- 
ure [12] shows the radial velocity of hot and cold gas ver- 
sus the radius at t=11.8 Myr for model 2D512xl28D (the 
data correspond to a snapshot shown in the right pan- 
els in Figure [7]). We reiterate that our model is quite 
simplified (e.g., no gas rotation) and the outer radius is 
relatively small (i.e., 200 pc). Therefore, our comparison 
is only illustrative. 

We find that the hot outflow originates at around 
10 pc and accelerates up to about v max m 200km/s, 
which is comparable to the escape velocity from 10 pc, 
v esc = 314 km/s. At larger distances, r = 100-200 pc, we 
see a signature of deceleration which is consistent with 
the observations of Seyferts outflows. We note that the 
geometry of the simulated flow is affected by our treat- 
ment of the boundaries of the computational domain, 



specifically along the pole and the equator we use reflec- 
tion boundary conditions (see §[3]). 

The scatter plot also indicates that the cold clouds ap- 
pear at about 20-80 pc. Their maximum velocity is about 
Vmax = 100 km/s, which is smaller than the velocity of 
the hot outflow. The plot does not show a clear indica- 
tion of a linear acceleration of the outflowing cold gas. 
However, it is possible that the cold clouds, seen in this 
snapshot, will continue to be dragged by the hot outflow 
and eventually will reach higher velocities. 

We also measured the column density of the hot and 
cold gas for the same, representative, snapshot at t— 11.8 
Myrs. The typical column densities vary with the ob- 
server inclination, Njj = 5 x 10 22 — 10 24 cm 2 for gas 
T > 10 5 K, and N H = 10 20 - 10 23 cm 2 for gas with 
T < 10 5 K. This is roughly consistent with column den- 
sities estimated from observations of AGN (e.g. for NGC 
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Fig. 8. — Evolution of the volume filling factor f vo i in model 
2D5f 2xl28D (black, solid line) and 2D256x64D (dashed, blue line). 
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Fig. 9. — Mass accretion rate in 2-D models with initially seeded 
random perturbations. Various colors code the M in models cal- 
culated with various grid resolutions: 2D256x64 (blue), 512x128 
(red), 1024x256 (green). Results are sensitive to the resolution 
same as in 1-D models. 



1068 N H = 10 19 - 10 21 cm 2 , iDas et all I20071 and refer- 
ences therein). 

Our results are similar i n many resp e cts, to the pre- 
vious findings presented in iBarai et al.l (|2012f) , i.e. the 
accretion evolution depends on fx luminosity; we also 
observe clouds, filaments and outflow. The outflow ap- 
pears at fx = 0.015 which is consistent with fx = 0.02 
found bv lBarai et al.l (|2012l) . Here, we are able to calcu- 
late models for about 10 times longer in comparison to 
3-D SPH models. We confirm the previous results that 
the cold phase of accretion rate can be only a few times 
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Fig. 10.— Model 2D512xl28D: In-, out- and total mass flow rate 
as a function of radius for three time moments shown in Figure 171 
(green, blue, and black correspond to t=3, 6 and 11.8 Myr). The 
solid lines mark the inflow rates while the dashed line - outflow 
rate. The dotted line is the total mass flow. 

larger in comparison to the hot one. 

Si milar models has been investigated in the past by 
e.g. IKrolik &; London! (|1983| ). Our work is on one hand 
a simplified and on the other hand an extended version of 
these previous works. The key extention here is that our 
new results cover the non-linear phase of the evolution. 
There are two new conclusions added by our analysis to 
the previous investigations. The 2-D models with out- 
flows are possibly governed by other than TI instabilities 
mainly convection. Another non-linear effects found in 
our 1-D and 2-D models is that the fragmentation of the 
flow makes it optically thick for photoionization. Further 
investigation of shadowing effects is required. 

Some sub-resolution models of AGN feedback in galaxy 
formation (iDi Matteo et al.l 120081 : iDubois et afl 120101 : 
ILusso fc Ciottill201lD assume that BH accretion is domi- 
nated by an unresolved cold phase, in order to boost up 
the accretion rate obtained in simulations. Our results 
indicate that the cold phase accretion is unlikely domi- 
nant as even in well-developed and well- resolved multi- 
phase cases, the accretion is typically dominated by a 
hot phase. However, we note that the cold phase of our 
solution might be an upper branch of some more compli- 
cated multi-phase medium (i.e., a mixture of molecular, 
atomic and dusty gas). 

This work was intentionally focused on a very limited 
number of processes and effects. Its results suggest that 
the future work should include more self-consistent ap- 
proach not only with shadowing effects but also with the 
radiation force. Our next step would be to investigate 
the non-axisymmetric effects via fully 3-D simulations. 
The latter is challenging and one may not be able to see 
very fine details of the gas dynamics as in 2-D models 
due to resolution effects. 
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Fig . 11. — Fraction of central illuminating source radiative energy 
intrinsically absorbed (upper panel) and emitted (bottom panel) 
by gas per second as a function of time in models 2D512xl28 (red, 
solid line) and 2D256x64D (blue, dotted line). 
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Fig. 12.— Scatter plot of radial velocity of hot (T > 10 s K, 
smaller red symbols) and cold (T < 10 5 K, larger blue symbols) 
phase of the flow in model 2D512xl28D at t=11.8 Myr (model 
shown in right panels in Figure [7)1 ■ 
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APPENDIX 



GROWTH RATE OF A CONDENSATION MODE IN A UNIFORM MEDIUM -CODE TESTS 

iFieldl (|1965() formulated a linear stability analysis of a gas in thermal and dynamical equilibrium. Here, we briefly 
recall his most important, for our analysis, equations. We disregard the thermal conduction effects. The dispersion 
relation derived from linearized local fluid equations with heating/cooling described by C function and perturbed by 
a periodic, small amplitude wave given by exp(nt + ikx), is: 



n 3 + N v n 2 + k 2 c 2 s n + N p k 2 c 2 s = 
where k is the perturbation wave number [k — 2tt/X) and functions N p and N v are defined as 



and 
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with c p and c v being the specific heats under constant pressure and constant volume conditions, respectively, and 
T is the gas temperature. Vertical line means that the derivative is taken under constant thermodynamical variable 
condition. Dispersion Equation IA1I has three roots. In a short wavelength regime (A <C 2ttN p /c s ), two, complex roots 
correspond to two conjunct nearly adiabatic sound waves and third, real one is an isobaric condensation mode (the 
gas density and temperature change in anti-phase so that the pressure remains constant). The sign of the real part 
of the root gives the stability criterion. The sound wave will grow if dC/dT\ s < (known as Parker's criterion, 
lParkerlfl953f ). The condensation mode will grow if dC/8T\ p < (Field's criterion). In a short wavelength limit, the 
growth rates asymptote to n = —0.5(N V — N p ) (for sound waves) and n — —N p (for condensation modes). Isochoric 
modes (n — ¥ —N v ) and effective acoustic waves are eigen modes of long wavelengths perturbations. The perturbation 
growth/damp time s cale is tti = 1/n. 

We use the above [Field! ([1965) theory to show that our numerical scheme for solving the modified energy conser- 
vation equation (Equation [3]) together with two other fluid dynamics equations is accurate. The test calculations are 
carried out in 1-D Cartesian coordinates within x € (0, L) range where L is the size of the computational domain in 
dimensionless units. The boundary conditions for all variables are periodic. In an unperturbed state, the gas density 
(po = 1) and internal energy density (eg = 1) are constant in the entire computational domain. The velocity of gas is 
set to zero. We assume that the gas is heated by an external source of radiation and cools due to free-free transitions. 
The test cooling function is simple: 

C = CpT 1/2 - H (A4) 

The normalization constants H (for heating) and C (for cooling) are set so that in the unperturbed state the gas is in 
radiative equilibrium i.e. £(po,eo) — 0. In this test the functions N p and N v have explicit, analytical forms 
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The numerical values of limiting growth/damp rates are N p = —0.04, and N v — 0.067, while the speed of sound is: 
c 2 = 1.11 (7 = 5/3). The domain sound crossing time is much shorter than the perturbation growth time scale which 
allows to keep the constant pressure. 

Our numerical scheme implemented into ZEUS-MP code correctly reproduces the expected growth rates of small 
amplitude perturbation of the uniform medium. The perturbation is an eigen mode of TI, and its properties depend on 
the assumed A. Eigen modes are realized by first applying a cosine pertu rbation to t he gas density p = po + Apg cos(kx) 
and calculating profiles of e and v from e .g. Equations 11 and 14 in (|Fieldl RL965T ) . for a given k and corresponding 
theoretical value of n (given by Equation lAlj) . Next we measure how fast the perturbation grows while it is in the 
linear regime. Figure [T3l (le ft panel), shows the analytical solution of the theoretical dispersion relation n(A) (solid line, 
third root of Equation IA1[) . and the numerical growth rates calculated with ZEUS-MP (points). For very short A's 
the eigen mode of this root is converging to the isobaric condensation mode and grows at n — —N p rate, as expected. 
The long A modes grow slower in comparison to the very short A condensations, as predicted by theory. For relatively 



14 



large A, the third root changes i nto an effe ctive acoustic wave, it becomes complex with the real part negative meaning 
that the waves are damped (see lShulll992l , Equation 41 in the Problem Set No 3). 
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Fig. 13. — Left panel: Analytical (solid line; third root of Equation lAlj! and numerical (points; calculated with ZEUS-MP) linear growth 
rates of an eigen mode with wavelength A. The initial amplitude of a perturbation is A = 10~ 3 . The numerical solution uses N x = 64 
grid points. Right panel: The condensation mode (density profile) with A = 0.1L shown at t=10 T (where T here is a sound crossing time 
over the computational domain L = 1) when it is resolved by 4, 8 and 16 points. The dotted line shows the very high resolution of 1600 
points for which the condensation mode grows at the expected theoretical rate. Here the growth rate is as expected when resolved by at 
least with 16 points. For lower resolutions the condensations grow slower than one would expect. 

In the second test, we measure the growth rate of a condensation mode that has a finite size (i.e. smaller than the 
domain length). We are interested in how many numerical grid points is required to resolve the correct n. We set 
A = 0.1 while L = 1. Figure [T3] (right panel) shows the same time snapshots of the growing condensation mode density, 
calculated with various numerical resolutions. Models with lower resolution evolve slower. When A resolved with 16 
points it starts converging to the right solution. We conclude that about 20 or more grid points per A is required to 
resolve the isobaric condensation. 



